home
***
CD-ROM
|
disk
|
FTP
|
other
***
search
/
Language/OS - Multiplatform Resource Library
/
LANGUAGE OS.iso
/
cpp_libs
/
rwvector.lha
/
RWVector2.1
/
src
/
dfft.cc
< prev
next >
Wrap
C/C++ Source or Header
|
1989-08-18
|
5KB
|
167 lines
/*
* Definitions for the double precision FFT server.
*
* Copyright (C) 1988, 1989.
*
* Dr. Thomas Keffer
* Rogue Wave Associates
* P.O. Box 85341
* Seattle WA 98145-1341
*
* Permission to use, copy, modify, and distribute this
* software and its documentation for any purpose and
* without fee is hereby granted, provided that the
* above copyright notice appear in all copies and that
* both that copyright notice and this permission notice
* appear in supporting documentation.
*
* This software is provided "as is" without any
* expressed or implied warranty.
*
*
* @(#)dfft.cc 2.1 8/18/89
*/
#include "rw/DoubleFFT.h"
#include "rw/mathpack.h"
DoubleFFTServer::DoubleFFTServer()
{
server_N=0;
}
DoubleFFTServer::DoubleFFTServer(unsigned N)
{
setOrder(N);
}
DoubleFFTServer::DoubleFFTServer(const DoubleFFTServer& s)
: (s), roots_of_1(s.roots_of_1), conjroots_of_1(s.conjroots_of_1)
{
server_N = s.server_N;
}
void
DoubleFFTServer::operator=(const DoubleFFTServer& s)
{
DComplexFFTServer::operator=(s);
server_N = s.server_N;
roots_of_1.reference(s.roots_of_1);
conjroots_of_1.reference(s.conjroots_of_1);
}
void
DoubleFFTServer::setOrder(unsigned N)
{
server_N = N;
DComplexFFTServer::setOrder(N);
roots_of_1.reference(rootsOfOne(2*N, N/2+1));
conjroots_of_1.reference(conj(roots_of_1));
}
/*
Performs DFT of a real sequence. Must have an even number of points.
The forward fourier transform of a real sequence is a complex
conjugate even sequence. If the original sequence had 2N points, the
result will have N + 1 complex points, for a total of 2N+2 points.
The extra two points are the imaginary parts of C(0) and C(N). Both
are always zero. Reference: Cooley, et al. (1970) J. Sound Vib.
(12) 315--337. This is algorithm #4. See the class description
header for more information.
*/
DComplexVec
DoubleFFTServer::fourier(const DoubleVec& v)
{
unsigned Nstore = v.length();
unsigned N = Nstore/2;
unsigned Nhalf = N/2;
unsigned Nhalfp1 = N/2+1;
fortran_int Nint = N;
// Various things that could go wrong:
checkEven(Nstore);
if(!Nstore) return DComplexVec();
if(N != server_N) setOrder(N);
#if COMPLEX_PACKS
// Bit of a cludge, in the interests of efficiency:
DComplexVec A((DComplex*)v.data(), N);
DCfftf_(&Nint, A.data(), the_weights.data());
#else
DComplexVec tempcomplex(v.slice(0,N,2), v.slice(1,N,2));
DComplexVec A = DComplexFFTServer::fourier(tempcomplex);
#endif
DComplexVec results(N+1);
DComplexVec A1(Nhalfp1);
DComplexVec A2(Nhalfp1);
A1.slice(1,Nhalf,1) = conj(A.slice(N-1,Nhalf,-1)) + A.slice(1,Nhalf,1);
A2.slice(1,Nhalf,1) = DComplexVec(imag(A.slice(N-1,Nhalf,-1))+imag(A.slice(1,Nhalf,1)),
real(A.slice(N-1,Nhalf,-1))-real(A.slice(1,Nhalf,1)));
A1(0) = DComplex(2.0*real(A(0)),0);
A2(0) = DComplex(2.0*imag(A(0)),0);
DComplexVec temp = A2 * conjroots_of_1;
results.slice(0,Nhalfp1, 1) = DComplex(0.5) * (A1 + temp);
results.slice(N,Nhalfp1,-1) = DComplex(0.5) * conj(A1 - temp);
return results;
}
/*
Take the IDFT of a complex conjugate even sequence. V is assumed to
be the lower half of the complex conjugate even sequence. See the
comments in the class header about the format of even sequences. The
results is a real periodic sequence. Reference: Cooley, et al.
(1970) J. Sound Vib. (12) 315--337. This is algorithm #5.
*/
DoubleVec
DoubleFFTServer::ifourier(const DComplexVec& v)
{
unsigned N = v.length()-1;
unsigned Ntot = 2*N;
unsigned Nhalf = N/2;
unsigned Nhalfp1 = N/2+1;
if(N != server_N) setOrder(N);
DComplexVec lowslice = v.slice(0,Nhalfp1,1);
DComplexVec hislice = conj(v.slice(N,Nhalfp1,-1));
DComplexVec A1 = lowslice + hislice;
DComplexVec A2 = (lowslice - hislice) * roots_of_1;
DComplexVec A(N);
A.slice(0,Nhalfp1,1) = DComplexVec(real(A1)-imag(A2), real(A2)+imag(A1));
A.slice(N-1,Nhalf,-1) = DComplexVec(real(A1.slice(1,Nhalf,1))+imag(A2.slice(1,Nhalf,1)),
real(A2.slice(1,Nhalf,1))-imag(A1.slice(1,Nhalf,1)));
DComplexVec X = DComplexFFTServer::ifourier(A);
#if COMPLEX_PACKS
DoubleVec results((double*)X.data(), Ntot);
#else
DoubleVec results(Ntot);
results.slice(0,N,2) = real(X);
results.slice(1,N,2) = imag(X);
#endif
return results;
}
void
DoubleFFTServer::checkEven(int l)
{
if( l%2 ){
char msg[60];
sprintf(msg, "Length (%d) of a real series must be even.", l);
RWnote("DoubleFFTServer::[i]fourier()", msg);
RWerror(DEFAULT);
}
}